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Abstract 

We study the problem of computing energy density in one-dimensional quantum systems. We 
show that the ground-state energy per site or per bond can be computed in time (i) independent 
of the system size and subexponential in the desired precision if the ground state satisfies area 
laws for the Renyi entanglement entropy (this is the first rigorous formulation of the folklore 
that area laws imply efficient matrix-product-state algorithms); (ii) independent of the system 
size and polynomial in the desired precision if the system is gapped. As a by-product, we prove 
that in the presence of area laws (or even an energy gap) the ground state can be approximated 
by a positive semidefinite matrix product operator of bond dimension independent of the system 
size and subpolynomial in the desired precision of local properties. 


1 Introduction 

The quantum PCP conjecture [2] is an assertion about the computational complexity of the ground- 
state energy per bond on general interaction graphs [6]. It is one of the most prominent open prob¬ 
lems for theoretical computer scientists in the emerging field of quantum Hamiltonian complexity 
|14j . Furthermore, computing energy density on regular lattices is a fundamental problem in con¬ 
densed matter physics. Here we present a comprehensive study of this problem in one-dimensional 
(ID) quantum systems. 

Consider a Hamiltonian H = ^ chain of n spins (qudits), where Hi acts on the 

spins i and f -|- 1 (nearest-neighbor interaction). Without loss of generality, assume \\Hi\\ < 1 
and X{Hi) = 0, where A(- • •) denotes the ground-state energy of a Hamiltonian. We would like 
to compute the energy density \{H)/n to precision 5, where 5 = 0(1) is an (arbitrarily small) 
constant. We begin with a general result. 

Lemma 1. A(i7)/n±(5 can he computed in time with probability at least 0.999. 

Proof. We divide the system into n6/2 blocks, each of which consists of 2/5 spins (we assume 
both n6/2 and 2/6 are integers for convenience). In particular, define H' = H[ with H[ = 

H 2 i/s+j acting on the spins 2i/6 + 1, 2i/6 -|- 2,... , 2(z -|- l)/5. Let ii,i 2 , ■ ■ ■ ,im with m = 


1 


0(1/5^) be independent uniform random variables on {1, 2,... , re(i/2}. The Chernoff bound implies 


\H - H'W = 


n6l2 


2=1 


< n5/2 ^ 0 < \{H)/n - \{H')/n < 5/2 


Pr 




k=l 


> <5 < Pr 




k=l 


> 5/2 


< exp(—n(m(i^)) < 0.001. 


( 1 ) 


As each can be computed by exact diagonalization, and the running time of the algorithm 

is mexp(0(l/(5)) = exp(0(l/(5)). □ 

Remark. This algorithm is an efficient polynomial-time randomized approximation scheme (EPRAS) 
but not a fully polynomial-time randomized approximation scheme (FPRAS). Here, a randomized 
f{l/5) poly(n)-time algorithm is an EPRAS (/ can be any function), and a randomized poly(n/(5)- 
time algorithm is an FPRAS. As computing \{H) to precision 5 = 1/ poly(n) is QMA-complete 
[3 SI, an FPRAS in general ID quantum systems is believed to be impossible. 

The algorithm of Lemma[l]is inefficient both in theory and in practice because its running time 
grows exponentially with 1/5. In this paper, we improve the scaling in 1/5 in the presence of area 
laws for entanglement or even an energy gap. 


2 Computing energy density in the presence of area laws 

A folklore in condensed matter physics is that in ID a (ground) state can be efficiently represented 
by a matrix product state (MPS) if it satisfies an area law for entanglement; then, it is likely 
that the MPS representation of the ground state can be found efficiently using the (heuristic) 
DMRG algorithm [19112f)j . In this section, we give the first rigorous formulation of this folklore. 
In particular, we show that the ground-state energy density X{H)/n can be computed to precision 
5 in time independent of n and sub exponential in 1/5 if the ground state satisfies area laws for 
the Renyi entanglement entropy (a function f{x) is subexponential if f{x) = o(exp(x'^)) for any 
constant c > 0). As a by-product, we prove that in the presence of area laws the ground state 
can be approximated by a positive semidefinite matrix product operator (MPO) of bond dimension 
independent of n and subpolynomial in the precision of local properties. 

We use the Renyi entanglement entropy because an area law for (or even logarithmic divergence 
of) the Renyi entanglement entropy Rq( 0 < a < 1) implies efficient MPS representations [18]. In 
contrast, an area law for the von Neumann entanglement entropy does not necessarily imply efficient 
MPS representations (see m for a counterexample), although the von Neumann entanglement 
entropy is the most popular measure of entanglement (for pure states) in quantum information and 
condensed matter theory. 

2.1 Preliminaries 

MPS is a data structure underlying DMRG-like algorithms. 

Definition 1 (Matrix product state (MPS) [151 [8]i. Let d = 0(1) be the local dimension of each 
spin, and {|jA:)}^j.=i be the computational basis of the Hilbert space of the spin k. Suppose 
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with Di = Dn+i = 1 is a sequence of positive integers. An MPS with open boundary conditions 
takes the form 

W= E (2) 

jl j'2vdn = l 

where is a matrix of size x Dfc+i- Define D = max{Dfc}^^| as the bond dimension of the 
MPS 

Remark. Any state can be expressed exactly as an MPS of bond dimension D = exp(0(n)), and 
an MPS representation is efficient if D = poly(n). 

Any MPS can be transformed into the so-called canonical form m with the same bond dimen¬ 
sion such that 

^ = 1, ^ = Al^+H, (3) 

j=i i=i 

where A^^l = diag{(A|^^)^, ...} with aJ^^ > A^^^ > ••• > 0 the Schmidt coefficients across 

the cut k — l\k in nonascending order. In the remainder of this section, we will use the canonical 
condition ([3l) extensively. 

MPO is the operator analog of MPS. 

Definition 2 (Matrix product operator (MPO)). Let be a basis of the space of operators 

on the spin k. Suppose with Di = Dn+i = 1 is a sequence of positive integers. An MPO 

with open bonndary conditions takes the form 

ild2,...Jn = l 

where is a matrix of size Dk x Dk+i- Similarly, define D = max{iAfc}^^| as the bond dimension 
of the MPO K. 

Remark. Any operator can be expressed exactly as an MPO of bond dimension D = exp(0(n)). 

Definition 3 (Entanglement entropy). The Renyi entanglement entropy Ra{0 < a < 1) of a 
bipartite (pnre) qnantum state pab is defined as 

Ra{pA) = {'^-p% (5) 

where pA = tr^ PAB is the reduced density matrix. The von Neumann entanglement entropy is 
defined as 

5'(m) =-tr(pAlog/OA) = lim Ra{pA)- (6) 

Remark. As Ra is a monotonically decreasing function of a, an area law for R^ implies that for 
Rp if a < (3. 

Fixing a cut, the truncation error can be upper bounded as follows. 

Lemma 2 ([T8]L Let Ai > A 2 > • • • > 0 with X'j = 1 be the Schmidt coefficients across the cut 
in nonascending order. Then, 

£d-= Y < exp((l - a)(R„ - logD)/Q:). (7) 

j>D+l 
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Minimizing the energy over MPS of constant bond dimension can be done efficiently by dynamic 
programming (or the transfer matrix method in statistical mechanics). 

Lemma 3 l[T6l[T]l. Let AiVS £) he the set of all normalized MPS of bond dimension D, and H be a 
local Hamiltonian on a chain of n spins. Then can he eomputed to precision 

1/poly(n) in time 


2.2 Approximation of local properties 

As errors may accumulate while truncating each bond, it is unlikely that the ground-state wave 
function can be approximated up to a small constant error (measured by fidelity) by an MPS of 
bond dimension independent of the system size, even if there is an energy gap. In this subsection, 
we prove that in the presence of area laws the ground state can be approximated by a positive 
semidefinite MPO of bond dimension independent of the system size and subpolynomial in the 
precision of local properties (a function f{x) is subpolynomial if f{x) = o{x'^) for any constant 
c > 0, i.e., it grows slower than any power of x). 

Lemma 4. There exists a positive semidefinite MPO pd of bond dimension sueh that 


{fi\d\fi) -tr {^pdO^ 


< 0{y/£D), Sd 


max 

k 


E 

j>D+l 



( 8 ) 


for any local operator O with bounded norm, where are the Schmidt coefficients of \ip) 

across the cut k — l\k in nonascending order. 

Proof. We express {fi) exactly as an MPS ([2|) of bond dimension exp(0(n)) in the canonical form 
([3]) . and define an MPO pn oi bond dimension as 


PD= • • • in){jlj2 • • • Jnl- (9) 

In the bulk, is a matrix of size x D^: at the boundaries, ■ and bY ■ are row 

’ ’ ’ *ljl ^n,3n 

and column vectors of length respectively. Without loss of generality, we will only work in 

ffcl 

the bulk in the remainder of this proof. The rows of S] ■ are labeled by two indices a, a for 
a, a' = 1,2,... ,D, and the columns of B^- are labeled by j5,j5' for j3,j3' = 1,2,... ,D. Let 

:= (A^^^)^ and S := {1,D + 1,D + 2 ,...}. Define xY = if s E 5 and xY = 1 

ffcl 

otherwise. The matrix elements of j are given by 


Bljia = a' = 1] (3 = fi' = 1) = ^ Af\s, t)Af^*{s, t) (x^^ , 

s,teS 

Bl^jia = a' = l;fi + fi'>3) = Af'(s,/3)Af'*(s,/?') (^x^f 

seS 

Bljia + a' = fi' = 1) = ^Af'(a,t)A^^'*(a', 

teS 

B^ka + a' > 3;/3 + /3' > 3) = Af'(a,/I)Af'*(«',/3')xL"'x£', 


( 10 ) 
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where (e.g.) t) is the element in the sth row and tth column of Dehne = X^teS 

if s G S' and = A^^^ otherwise. The matrices satisfy a pair of conditions similar to (j3]) : 

B|5(a,«';/?,/3')V = -5-', E (H) 

/3,/5',i a,a',i 


Indeed, 


^i?g](a = a' = l;/3,/3')V= E ^ 


/3,/3',i 


sGS,t,i 


ses 


E + a' > 3;/3,/3') <5/3/3' = Af\a,t)Af'^* {a',t)x^a^x^a} = = daa', (13) 


/3,/3',i 


t.i 


E Bf}{a,a'-I3 = /3' = = E {s,t) 


a,a' ,i 


t£S,s,i 


Y: Af'(s,t)Af%,t)AP =Af*", 


(14) 


t(^S,s,i 


Bf;}ia,a'-f3 + f5' > 3)d„„,AM = Afl(.,/3)Afl*(s,/3') (xW)Em 


a,a',2 


= Y4\s,^)AP*is,l3')A^^^ =6p^:A 


,[fc] 


[fc+i] 


(15) 


Assume without loss of generality that O is a 2-local operator on the spins k — 1 and k. Let 
{4^^)}f=i be the computational basis of the Hilbert space of the spin k. The matrix elements of 
the reduced density matrix of on the support of O are 


where 


Let 




^[fc-l] 7[A;] 




s',t'E5,s,£>l 


[k] 

p 


q+/3>3,s,£ 

— *0 *0 > 


s,t,a,IS 


(16) 


Y Af ^\s,a)Af ^^*{s,f5)A^^ ^^Af\a,t)AY{f3,t)x]^^x^^\ (17) 

E Af-^\sA')Af-^^HsA')Af-^^A^\s'A)AYis'A) • (18) 


s' ,t'GS,s'+t'>3,s,t 




(19) 
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such that 


*1 ~ 


-j _ 2 

riti'* ^ ^ ^ \ ^ ri-i' , rid' 






E l ^ t'^ 4' r'^4'* 

I'^l — ^s,t,s',t'^s,t,s',t' 

' JJ' ,jj' 

s' ^t'£S^s'+t'> 3 ,s^t,i,i' 

O r/2 A [fc] A 

^ a!!=’(xW) = ^s' ^^^2 ^ Al^]<2dhD. (20) 


J l^l J l^l =f= 

^s,t,s',t'-‘^s,t,s',t' 


E Af'hl';') = 

s' ,t' £S,s'+t'>3 


YlreS 


s>D+l 


Let P = diag{l, 1,..., 1, 0,0,...} with tr P = D and Pi = diag{l, 0,0,...} < P. The matrix 
elements of the reduced density matrix Qf |•0)(•0| on the support of O are 


_ ^j,j'',iP _j_ ^jd'dP _j_ j,jd'\iP 


where 


^ tr - P)Af^Af^^ . 


Using the Cauchy-Schwarz inequality, 






X / ^ tr - P)Af^ 


< ^Jd^tT{PAi’^^P)d^tT{{l - P)AW(1 - P)) = d^ lY^ Ai^^ Y I Y1 ^ 


s<D s>D-\-l 


s>D+l 


< d^y/eB- 


Similarly, 


r^’-^ ’*’* < \Jd? tr A^d^ tr((l — P)A[^1(1 — P)) = / Y^ — d^V^- (^4) 

V s>D+l 
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Using the Cauchy-Schwarz inequality, 


(s, a)Af-^^*{s, /3 )Af {a, t)A^^* (/3, t) 

a,^<D,s,t 


E 


1^0 


’0 


1 _ 


< Y. 4""'^ a) Af(s, /3) W (a, t)Af* (/?, t) 

a,l3<D,s,t,i,i',j,j' 

< Y Af-^\s,a)Af-^^*{s,/3)Af-^^Af\a,t)AYW,t) (l - + 1 - 

a,/3<D,s,t,i,i' ,j,j' 

J- / / ^ ^ ^ 3 3 * '^3 


< (l - xf') cf /Afl a'^ 1 + /Afl Ai^l < 2dVAfl (l - xf' 

\ V V / 

= 2rf=Ari (^1 - Afl /YfY “ ^ "" \/W + "°) 

< 2d2e^ j (2aP +eD/AS^') < 

Finally, ([8j) follows from 


(25) 






+ 


+ 


I n *n 


+ 


< d‘^y/££) + d‘^y/££) + d‘^y/2££) + 2d3££) — 0 (\/ £d') 

= O(v^). 




as d = 0(1) is an absolute constant. 


(26) 

□ 


Theorem 1. Suppose jV') satisfies an area law for the Renyi entanglement entropy Ra (across any 
cut). Then there exists an MPO p of bond dimension exp(2PQ,)0(l/d)^"/(^“") such that 


{fi\d\fi) - tr {^p6^ 


< 5 


(27) 


for any local operator O with bounded norm. 

Proof. This is an immediate consequence of Lemmas [2l[31 In particular, 

S = 0(y/^) < 0(exp((l - a)(Ra - log D)/(2a))) => = exp(2P„)0(l/d)^"/(i-“). (28) 


□ 


Remark. If Ifi) satisfies area laws for the Renyi entanglement entropy R\/a, then the bond dimension 
is independent of n and subpolynomial in 1/5. 
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2.3 Algorithm and analysis 

Theorem 2. Suppose the ground state |'0) of a loeal Hamiltonian H on a chain of n spins satisfies 
an area law for the Renyi entanglement entropy (across any cut). Then there exists an MRS 
\4>) of bond dimension D = exp(i?Q,)0(l/(i)^“/d-«) such that 


{(t)\H\(f) < \{H) + n5. (29) 

Proof. Recall that pD dH]) is an MPO of bond dimension such that for any local operator O 
with bounded norm, 

tr {j)dO^ < (V’|0|V^) + 5 ^ iT{pDH) < (V’l-fflV') + n5. (30) 

The matrices defined in (fTOl) can be expressed as a sum of many terms: B^f] = where 

c'fjpq is a matrix of size x (except at the boundaries). In particular, the matrix elements 
ffcl 

of Cl j p^^ are given by 


^£Uo,5>i(«> /3 = /?' = 1) = df'(a, D + g)df *(«', D + (32) 

^£Ui,.=o(« = «' = = Af\D+p,l3)A^f*{D+p,l3') {JS+p)\ (33) 

^£Ui, 5 >i(« = o' = 1; /3 = /?' = 1) = Af\D+p, D + q)Af*{D+p, D + q) (xg^ ,(34) 

^5Uo,,>i(«>/3 + /?' > 3) = c£],>i,,=o(« + «' > 3; /3, /3') 


= ^£Ui,.>i(«>/3 + /3' > 3) = C‘';‘p>i,,>i(a + a' > 3; /3, /3') = 0. 


[fc] 


(35) 


Thus, the MPO pD can be expressed as a sum of many MPO of bond dimension D^: 


PD 


E 


^^,qi,P2:q2,...,Pn,qn 


Pl,qi,P2,q2,---,Pn,qn 


(36) 


where 


^i,qi,P2,g2,...,Pn,gn 


h ,il,*2,i2,-..,*n ,in = l 


/^[ 1 ] ^[ 2 ] _ _ _ 
V^*l,il,Pl,gA*2,j2,P2,g2 ' ' ' ^in,jn,Pn,qn 


hi2---in){jlj2---jn\- 


(37) 

We observe that each p^’'Ji:P 2 ,g 2 ,...,Pn,qn -g density matrix of an unnormalized (pure) MPS of 
bond dimension D. Hence (j30p implies that one of them (still denoted by p^’'Ji’P 2 >'? 2 , -,Pn,gn^ satisfies 


^.^^pP^,qi,P2,q2,..;Pn,qn ^^,qi,P2,q2,...,Pn,qn 

and we complete the proof by normalizing p^’'Ji:P 2 ,g 2 ,...,Pn,qn^ 
Let 0(x) := 0(x poly log x) hide a polylogarithmic factor. 


< (V'|77|V’) + n6, 


(38) 

□ 


Theorem 3. Suppose the ground state \^p) of a loeal Hamiltonian H on a chain of n spins satisfies 
an area law for the Renyi entanglement entropy Ra (across any cut). Then \{H)/n ± 5 can he 
computed in time exp(exp(2i?Q,)0(l/(5)^"/d-«)) y^ah probability at least 0.999. 


Proof. We divide the system into n6/2 blocks, each of which consists of 2/6 spins. In particular, 
define H' = with ^ H 2 i/s+j acting on the spins 2i/6+l, 2i/6+2 ,... , 2{i + l)/5. 

Theorem [2] implies an MPS |0) of bond dimension D = exp(i?Q,)0(l/5)^“/(^““) such that 

<\{H)+n6/2. (39) 


The Schmidt decomposition of |</>) across the cut 2i/5\2i/6 + 1 implies 
D D 

\<t>) = Y,HL,)\R,) ^ = ^ s. t. < m'M)- (40) 

i=i i=i 

As \Rki) is an MPS of bond dimension D, the Schmidt decomposition across the cut 2(i + l)/<5|2(i + 
l)/<5 + 1 implies 


D 


D 


\Rk,) = ^ {Rk^Rk/) = 

i=i i=i 

^s.t. (L'y\H[ < {Rk^Rk:} < {m'M) 

n6/2 n5/2 n5/2 

^ \{H) - n6/2 < X{H') = ^ A(i/') < ^ (l',, i/' L/,) < 


i=l 


i=l 


i=l 


= (())|i/V) < < X{H)+n6/2. 


(41) 


As \L',,) is an MPS of bond dimension D, Lemma [3] implies that a number between X{H[) and 
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{L'y\H'/\L'y) can be computed in time {1/6)^^^ \ Finally, a randomized algorithm with running 
time 0{l/6'^){l/6)^^^^^ = exp{ex.p{2Ra)0{l/6)'^°‘^^^~°‘^) follows from the proof of LemmalU □ 

Remark. If {fj) satisfies area laws for the Renyi entanglement entropy R\/a, then the running time 
is independent of n and sub exponential in \/6. 


3 Computing energy density in the presence of an energy gap 

3.1 Introduction 

As an energy gap implies area laws for the Renyi entanglement entropy Rvo m, all results in 
Section [2] apply. In particular, to approximate local properties the bond dimension of the MPO is 
subpolynomial in l/d (Theorem [1]) , and to compute the ground-state energy density the running 
time of the algorithm is sub exponential in l/d (Theorem [3]). Indeed, the scaling in l/d can be 
obtained precisely. 

Let I To) and e be the ground state and the energy gap of H, respectively. Since we are 
interested in the regime l/d is large, we assume 1/6 S> 1/e. Some e-dependent subpolynomial 
(e.g., 2^(^ i/^iog3/4(i/5))^ constant (e .g., 2^*^^^^'^) factors will appear below. If not dominant 
(e.g., accompanied with poly(l/5)), depending on the context they may be neglected or kept for 
simplicity or clarity, respectively. 

Let n(x) := D.{x/ poly log x) hide a polylogarithmic factor. As a by-product of the proof of the 
area law for entanglement, the truncation error can be upper bounded as follows. 
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Lemma 5 ([lO]). Let {A*} be the Schmidt coefficients o/|'I'o) across a cut. Then, 

^^20(iA+,-i/4iog3/4(,/,^))^ ^ (42) 

i>D+l 


Then, it is straightforward to obtain a pair of corollaries. 

Corollary 1. There exists a positive semidefinite MPO p of bond dimension 
such that 


(V'lOIV’) -tr 


< <5 


20(lA+e-l/4 log3/4(l/5)) 

(43) 


for any local operator O with bounded norm. 

Corollary 2. \{H)/ndcd can be computed in time 2^°^' '°® with probability at least 0.999. 


Our goal is to improve the scaling in 1/5. The main result of this section is 

Theorem 4 (FPRAS for computing ground-state energy density of ID gapped Hamiltonians). 
\{H)/n ± 0{5) can be computed in time (1/(5)*^^^^ with probability at least 0.999. 

Proof. Recall H = where 0 < Phi <1 acts on the spins i and i + 1. Following the proof of 

LemmalU we divide H the system into n5 blocks, each of which consists oil/5 spins. In particular, 
dehne H' = hT'- = Y^yJi ^ ^-^ting on the spins i/5 + l,i/6 + 2,..., {i + l)/5. 

Then, it suffices to prove Lemma [3 □ 


Let PL be the Hilbert space (C'^)®"' of the system (n qudits) and ^[ 44 ^] = (C'^)®^*’' be 

the Hilbert space of the spins with indices in the interval [ii,ir\- Let cl = hPj) and 

= KE]Zl+s+i Hj). Define 

Hp := H, - ci j Pf + t(l - Pf), (44) 

where P^* is the projection onto the subspace spanned by the eigenstates of Hj with 

energies at most cl + t-, similarly we define . Let 

P-* := Hj^ -I- Hij—s + Hi^—s+i + • • • + Phir+s + Hp ^ ir — ii P 2s P 2t P 1. (45) 


be a ID Hamiltonian which is local on P[jj_s_|_i 4 ^_|_s]. Let |Tq) and e' be the ground state and the 
energy gap of P-*, respectively. 

Lemma 6 ([lO])- For t > 0(log e“^), 

(a) 0 < X(H) - X(H^P -cl-crP 2-^W; 

(b) |('I"o|To)| > 1 - 2 -^W; 

(c) e' > e/ 10 . 


Definition 4 (support set). For any interval [ii,ir], S C Pl[^^k] with /c < P is a {k,s,b, Af> or Ay)- 
support set if there exists a state {i/j) £ PI (called a witness for S) such that 

(i) the reduced density matrix of |V^) on Pl[ii^k] is supported on span S'; 

(ii) |S| < s; 

(iv) all elements in S are MRS of bond dimension at most 6 ; 

(hi) |(V'|4 'q)| > 1 — Ay or < X{H-^) P Age' (depending on the context either Ag or Ay 

is used as the precision parameter). 
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Table 1: Evolution of the parameters in each iteration. The asterisks mark the parameter that is 


reduced at every step. 



k 

s 

b 

^.f 

Ae 

start 

k-1 

PlP3 

P 2 P 3 

n/a 


extension 

k 

dpips 

P 2 P 3 

n/a 

ce'i4 

cardinality reduction 

k 

Pi* 

dpip2pl 

1/1000 

1/1000 

bond truncation 

k 

Pi 

P2* 

1/20 

n/a 

error reduction 

k 

PlP3 

P 2 P 3 

n/a 



Lemma 7. A number e* between and + 3 can be computed in time 

Proof sketch. Set ii = i/5 + 1 and i^ = {i + ^)/5. Our algorithm is a “local” variant of and thus 
very similar to the polynomial-time algorithm for computing the ground-state wave function in 
ID gapped systems naiiniE]. In particular, it iteratively constructs a {k,pip^,p 2 P‘s, = ce'^^)- 
support set for A: = i/, z;-|-1,..., V, where Pi,P 2 ,Ps are (upper bounded by) ^-independent and e- 
dependent polynomials in 1 /5, and c is a sufficiently small absolute constant. After the last iteration, 
we solve the convex program = mintr(iTff^), whose variable u is a density matrix on spanS"*^. 
This convex program is of polynomial (in l/d) size because (i) spanSj^ is of polynomial dimension; 
(ii) any element in spanS'j^ is an MPS of polynomial bond dimension. Clearly, e* > X{Hl). By 
dehnition, there exists a state \'tpi) G H such that (i) + ce'’^ < X{H-^) -t- 1; (ii) 

its reduced density matrix on is supported on span5*^. Property (ii) implies e* < 

and (i) implies {i/ilH'-l'i/i) < A(ff') -|- 3. Indeed, let \4>i), {(/r) be the ground states of H/f = 

= Y!jX+i respectively. Then, 

+ {4>^\H'^4>i) + {4 >r\K\4>r) + ‘^> {(I>Mr\h^^\(/Mr) > A(^-') 

> - 1 > (z/ijliJ+lz/ij) -I- - 1 

> |0l) + + {4>r\H/^\4>r) - 1. (46) 

Each iteration consists of four steps: extension, cardinality reduction, bond truncation, 
and error reduction. Table [T] summarizes the evolution of the parameters s,b,Af or Ae in each 
iteration of our algorithm. 

The analysis below gives 

p3 = il/5f^^\ (47) 

The running time of the algorithm is a polynomial in Pi,P 2 ,P 3 so that Lemma [7] follows. The 
remainder of this section explains each step in detail. □ 

3.2 Preliminaries 

Lemma 8. (a) ('0|77-*|V’) < X{H-^) + rje' implies |('0|'kg)| > |(V'|^o)P > 1 — 

(b) |(V'|<?^>i)| > 1 - and \{'i/\4>2)\ > 1 - mjmply \{(l)i\4>2)\ > 1 - 2(r/i +m}; 

(c) |(V'|</’)| > 1 — ^ implies |(V’|0|V’) — (<?i|0|(/>)| < 2y^ for any operator O with ||0|| < 1. 

Proof. Just a few lines of algebra. □ 
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Definition 5 (truncation). Let |V^) = Schmidt decomposition of a state 

IV') G % across the cut i\i + 1, where the Schmidt coefficients are in nonascending order: Ai > 

A> • • • > 0. Define trunc^ = Yl^=i 

Lemma 9 (Eckart-Young theorem). The state {ip') = trunc^i |^)/|| trunc^i |^)|| satisfies {ip'\ip) > 

I(0|V')I for any state \(p) G H of Schmidt rank D (across the cut i\i + 1). 

Lemma 10 ([12]). Suppose \(fi) ^Ti is a state of Schmidt rank D (across the cut i\i + 1). 

|(truncz)/,,V’l<^)l > I(V'|0)IVr/> 0 ,V’G^. (48) 

Fix a cut i\i + 1 with ii <i <ir — 1 ra the definition (|45h of 

Lemma 11. Sets > 0{e~^)+0{e~^/‘^\o^^‘^{l/r])). Then, (V'I'^'q) > l—rj for\fii) = trunc^^ |'I'q)/|| truncs^ l'I'o)ll? 
where By^ = 20(Ve+e-^/Dog3/4(i/^)) _ 

Proof. This is a by-product of the proof of the area law for entanglement m- □ 

Lemma 12 ([H]). + pe' and rj < 1/10 imply (ip'\H-^\ip') < + 25y/rj 

for \Tp') = truncB^ \fi)/\\ truncg^ \ fi)\\. 

3.3 Algorithm and analysis 

In the {k — i;)th iteration the algorithm constructs a {k,pip 3 ,p 2 P 3 , Ag = ce'^^)-support set Sk from 
a {k — l,pip 3 ,p 2 P 3 , Ae = ce'^'^)-support set Sk-i returned in the previous iteration. Extension is 
trivial (see [H]), and constructs a {k,dpiP 3 ,p 2 P 3 , = ce'^^)-support set S^\ 

3.3.1 Cardinality reduction 

Let trjjj] p denote the partial trace over TL^ij] of a density matrix p on Pi. 

Definition 6 (boundary contraction). Let \'fi) = = J2f=i 

be the Schmidt decompositions of a state \'tp) G Pi across the cuts ii — l\ii and k\k + l, respectively. 

Let {|j)}^i be the computational basis of C^. Let Ur{iP) : Pi\^k+i,n\ be the isometry 

specified by UR{fi)\j) = \rf^^^) such that = Y^f=i G Piii^k] ® C^. The 

right boundary contraction Cr{'iP) is a density matrix on Pi[k,k] ® C^: 

CRifi) := (49) 

Similarly, let Uiffi) : —)• Pi\^l^i^_^ be the isometry specified by UL{'ip)\j) = such that 

^ ® The left boundary contraction is a density 

matrix on (g) 

C'l(V’) := t/^nV’)tr[i,+i,n](|V’)(V’l)b^L('0)- (50) 

Let Hl := Hf + ■= EU^^^Hr := + Hf. Define Hf = 

Hr - X{Hl), H'l^ := Hm - X{Hm),H'^ := Hr - X{Hr) so that X{Hf) = X{H'^) = X{H'^) = 0. 
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Lemma 13. Let p he a density matrix on <8> ® C^, and {ip) = = 

'^f=i be the Schmidt decompositions of a state \^p) G Ti across the cuts ii — l|i; 

and k\k + 1, respectively. The density matrix p' := UL{'ip)UR{pp)pUfp^{'ip)Uff^{ip) on TL has energy 

ti{p'H^^) < tv{pHM) + mHL + H,^-i + Hk + Hum 
+ P - Cl{\'iP))\\i I 1+ max ||i^ilOII 

y |Z>espan{|zWi)} 

+ l|trs^[i,,fc-i]/)-CR(|V’))||i I 1+ max ||i^ij|r')|| 1 , (51) 

y |r>espan{|rf+«)} j 

where p denotes the partial trace over TLi^^i^k] ® C^, and similarly for p. 

Proof. Write tic{p'H-^) = ix{p'HM) + tic{p'{HL + + ix{p'{Hk + Hr))- For the first term, 

ii{p'HM) = t,i{pHM) as UL{'tp),U r{'iP) are isometries. For the third term, 

tx{p{Hk + Hr)) - {'4)\{Hk + Hr)\\P) = tr[(/)' - \ii){'ip\){Hk + Hr)] 

= tr[(/9' - \'ip){'ip\){Hk + H'^)] = tr[tr[i_fc_i](/)' - \'ip){'ip\){Hk + H'^)] 

= tx[U])^{il!)tx[i^k-i]{p - \f^){f^\)UR{f>)U))^{il!){Hk + iL^)C/ij(V')] 

= tr[(trB®[ii,fc-i] P - CR{'ip))Uf^^{ip){Hk + H'^)Ur{'4))] 

< II p - CRm\i\\UR\mkUR{f^)+uf,\m'RURm\ 

<l|trs^[i,,fc_i]/)-CR(V^)||i 1+ iriax ||i^H|r)|| . 

y |r->espan{|ijf''+y)} j 

Similarly, for the second term, 

tx{p'{HL + Hg-i))-{'ip\{HL + Hi^-i)\'ip) < ||tr[ij+i^fc]^B/5-C'L(V')||i ( 1+ max 

y I0espan{|4**') 

Let Ni and W be, respectively, .^-nets with ^ = fl(e) for the trace norm over the space of left 
and right boundary contractions of bond dimension ~ 20 (iA) such that \Ni\ = |W| = 

It is straightforward to construct Ni and W m time poly|W| = poly |W| = 

220(1/') 


(52) 



□ 


Cardinality reduction convex program and bond truncation 

0. Let the variable p be a density matrix on C^iO'/'"'® (8)span5j.^^ ® C Q^w^e'^ ^ 

Q-^lOv^e'e _ 

1. For each Xi G Ni and each X^. G W; solve the convex program: 
min ixipHu); s. t. || P --^z||i < ?, II P --^r||i < ?, trp=l, p > 0. (54) 
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2. Let \ip) = Yliji \ ji)\V’ji,jr)\jr) be the eigenvector of the solution p with the largest eigenvalue. 

3. Let \(p') = ^ • ■ \ji)\^'A ,• )\jr) be the state obtained by truncating each bond (in whatever 
order) of \(p) to p 2 - 

4. Sj. consists of the MPS representations of all 

Let Pf be the projection onto the subspace (8i^[A:+i,n])“* spanned by the eigenvectors 

of H'^ + with eigenvalues at most t', and let Qt' be the projection onto the subspace spanned 
by the eigenvectors of with eigenvalues at most t'. 

Lemma 14 (truncation lemma). 

11(1 - PtOl^o)ll < 11(1 - Qi')l'I^o)ll < 100 • 2-*'/20. (55) 

Proof. The first inequality is obvious: Pt' > Qt' as = 0 and > 0. The second 

inequality was proved in m- □ 

Let t' = 0(log(l/e)) such that 100 • < ce'^^. 

Lemma 15. There exists a state \if) E ® span^i^^ (8) sueh that (i) its Schmidt 

ranks across the cuts ii — l\ii and k\k + l are upper bounded by (ii) (V’l 

250ci/®e'. 

Proof. Let |</>) be a witness for Since is a {k, dpiP3,P2P3, Ag = ce'^^)-support set. Lemma 
[81(a) implies |(())|'I'q)| > 1 — ce'^^. Lemma |8](c) implies 

+ RkM > + Hk)\%) - 

^ + H'j^ + < (ToKiL^ + H'j^ + Ll^)|4'o) + 4-\/^e'^ + ce'^"^. (56) 

Let \4>') = Pt/|(/>)/||Pt/|(/>)|| so that [(/)') E (^[i,*;-i] ®^[fc+i,„])-*'® spanC ®spanS’|,^^ ® 

'^^+1 n] construction. 

|(</>'|'L')| > \{m'\%)\ > |(</>|'L')| -|(</-|(l-PiOI'bo)l > l-ce'i"-||(l-Pt0l^'o)ll 

> 1 - - 100 • 2-‘'/20 > 1 - 2ce'^^ 

^ + HkM') < (T'o|(F,,_i + Ffe)|T'o) + 4V^e'" 

('I'ol(Lf)^ + + 6\/ce'’^ > + H'j^ + H'pf)\(f) 

= + H'^ + H'^)Pt'\4>) + mi - Pt'){Hf + H'j^ + - PuM) 

> mmL+ h'm++ t'ii(i - Ptmf im 

Clearly, X{Hf +Hh-i +H'm + Hk +H'j^) < \\m_,+Hk\\ < 2. Hence, H'm + KWo) < 2 

as > 0, and 

t' > + 6\/ce''^ 

^ {4>'\{Hl + P'm + P'b)W) ^ ('I'oK-f^L + P'm + Pr)\^o) + 6\/ce'^ 

^ + 6Vce'^ + 4Vce'^ < + lO^ce'^ (58) 

LemmafT^imolies that the state {(p”) := truncBjj,^^,6 l'?^^)/ll ^ 

^^+1 n] (1^® truncation is across the cut ii — l|i;) has energy {(p"\H-^\(p") < X{H-*) + IDDc^Pe'^. 
Once again, the state IV') := truncB^^^^j/^^,^ \^”)\\ ^ ® span5|,^^ ® 

'^fk+i n] (1^® truncation is across the cut k\k + 1) has energy {ip\H-*‘\ip) < X{H-^) + 250c^/®e'. □ 
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Lemma 16. S'® is a {k,pi,dpip 2 p‘^, Ag = 1/1000)-support set, where S® consists of the MPS 
representations of all \pji,jr)- 

Proof. Since Ni and Nr are ^-nets, there are elements Xi G Ni and Xr € Nr such that || — 
Cl('0) - Xi\\i < ^ and ||C'i?(V’) - ^r||i < f., respectively. Clearly, ii{pHM) < as 

UL{'if)UR{'i/)\f)){'ip\Uff^ is a feasible solution to the convex program (fMD . Let a = 

and set t = Q{e) such that 4^(1 + t') < eYdOOO. Lemma [T3] implies 

< tripHu) + + H,,-i + Hk + HR)\^|;) 

+ l|tr[i,+i,fc]^ijp-C' l(V')||i 1+ max \\Hf\l)\\ 

\ |0espan{|«'‘*')} J 

+ l|trB^[i,,fc-i]/5-C'R(V’)l|i 1+ max ||i?fl|r)|| 

y |r>espan||rj + J 

< {ijj\{HL +Hi^_l +Hm + Hk +Hr)\'4>)+ 12+ max \\Hf\l)\\+ max ||iSR|r)|| | 

< + 4^(1 + t') < A(44-*) + 250c^/^e' + e74000 < A(44-*) + e72000 (59) 


for sufficiently small constant c. We observe that 

(1) there exists at least an eigenstate of a with energy (with respect to 44-*) at most A(44-*) + 
e7l000; 

(2) there is at most one such eigenstate as Lemma [SKa) implies that such an eigenstate is close to 


4'; 


0/1 


(3) this eigenstate (denoted by |$)) has the largest eigenvalue due to Markov’s inequality in prob¬ 
ability theory; 


(4) |4>) = f7z,(7)C/i{(7)|<^) G ®'Hfk+i^n] is a witness for S'7'’ as a {k,pi,dpip 2 pl, Ae = 

l/1000)-support set with pi = ~ 


( 2 ) 


= 2 


20 ( 1/0 


□ 


3.3.2 Bond truncation 

Lemma 17. There exists a state lip) such that (i) |(?/)|'I/q)I > 1 — 1/1000; (ii) the Schmidt rank of 
17) across the cut j\j -|- 1 is upper bounded by r = 2*^i® ii < j < k — 1. 

Proof. We express IiI/q) exactly as an MPS ([2]) in the canonical form ([3|) with bond dimension 
exp(0(n)). Let P = diagjl, 1, ... , 1,0,0, ...} with trP = r. We define the unnormalized state |7*) 
as 


10) = E 


A'+f 

J1 J2 


.A^h-i]Mpj+i+i]p__PA 

3il /O + l 3k —1 


31,32,■■■,3’ 


.=1 


[fe —1] T} A[k] A [fc + 1] 


PAf'M 

3k Jfc+l 


Af%j2---3n). (60) 


Clearly, the Schmidt rank of |7*) across the cut j\j -|-1 is upper bounded by r for any ii < j < k — 1. 
Let E^'^\X) := ^ completely positive trace-preserving (CPTP) linear map. 

Then, 


(7'|4/'o) =trXfc+i, X,+i ■.= E^^\XjP), j = ii + l,ii + 2,...,k, Xy+i :=A[**+*1. (61) 


15 


As a CPTP linear map is nonexpansive with respect to the trace norm, 


||i?W(A)||i < ||A||i ^ IIA^+^l - A,-+1111 < ||aW - A,-Pill < ||A[^'1(1 - P)||i + ||A[^'1 - A,-Hi 

^ 1 - (V^'l^o) = tr(A['=+^l - Afc+i) < ||A['=+^1 - Afc+i||i < er{k - n) < er/5 < 1/1000, (62) 

where £r = (5/1000 for r = 2^^^^ <5 i)_ Similarly, 

(V>'|V’')=trA(,+i, A'+i :=P[^](PA'P), j = + 1, + 2,..., A:, A'+i := aI**+i] 

^ < ||A(,+i||i < • • • < ||A'+i||i = 1. (63) 

Hence |z/) := |V'0/IIIV')ll satisfies (V’I'I'o) > 999/1000. □ 

fs) 

Lemma 18. 5/ is a {k,pi,p 2 , Aj = l/20)-support set. 

Proof. Since 1$) is a witness for with energy |(‘h|P-*|$)| < A(P-*) + e'/lOOO. Lemma [SKa) 
implies |(<h|'I'Q)| > 999/1000. Lemma[8Kb) implies |(‘h|V’)| > 249/250. Let <h' be the state obtained 
by truncating the each bond j\j + 1 for j = ii,... ,k — 1 of ^ to p 2 = 1000r/(5. Then Lemma [TOl 
implies |(z/)|<l>')| > 199/200. Lemma[8Kb) implies |(4 'q|< 1>')| > 49/50. Hence |<h') = ULi'f’)UR{ip)\ip') 
is a witness for 5/ as a {k,pi,p 2 , Af = l/20)-support set. □ 


3.3.3 Error reduction 


Let p = ce'^^. Assume for the moment that we have an estimate eg of A(P-*) in the sense that 
i := |eg — A(P-*)| < t'!^/q < e'/2, where q = 41og(l/7y) + 24. Let 


A 



2 e '2 J 


y/27rq 



( e'^r^ 

V 2g 





dr. 


(64) 


Recall that 1$') is a witness for such that |(‘hQ|4>')| > 19/20. 

Lemma 19 ([H]). ((/5|P-*|(/)) < A(P-*) + r/e'/lOO for {(j)) = A|<1>')/||A|<1>')||. 


Let K be an MPO of bond dimension D across the cuts j\j + f for j = ii — l,ii,... ,k (its bond 
dimension across other cuts may be very large). It is straightforward to decompose A as a sum of 
p2 terms: 

D 

K= Y, ''4 ® . ( 65 ) 

jl,jr = l 

where are MPO on P[fe+i,n]) respectively, and is an MPO of bond dimension 

D on 

Lemma 20. There is an MPO Kt such that (i) ||Ar~exp(—zP-*r)|| < (ii) the bond dimension 
across the cuts j\j + 1 for ii — 1 < j < k is upper bounded by 2^^^^ j poly((5|') ; / Hi) ’s in the 

decomposition ([MP for Kt can he constructed in time f poly((5^'). 


Proof. We give a constructive proof so that (iii) will be apparent. Recall that H-^ is local on 
and has bounded norm on the left and the right of this region. Following [TB], for 
each ii — sl2 < j < ir + s(2 consider the unitary operator 




ir+S 


V{T) := exp -iH^^T - i ^ HiT exp -i ^ HiT - iH^^T exp(iP^*T). 


( 66 ) 


l=H-s 


l=j + l 
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It is straightforward to verify that V{T) = T exp i L(T)dTj is the time-evolution operator 
of the “Hamiltonian” L(t) = exp(—iII-*T)ffj exp(iff-^T), where T is the time-ordering operator. 
Define V'(T) = T exp i L'(T)dT'j with 

/ i+s'-i \ / i+s'-i \ 

L'(t) = exp j —i IIiT I ffjexp j i Hit | , s' < s/4. (67) 

\ i=j-s'+i j y i=j-s'+i j 

Since the “Hamiltonian” L'{t) acts on V'(T) also acts on Moreover, 

as a consequence of the Lieb-Robinson bound (see Ref. 24 in m for a simple direct proof) 

||L(r) - L'(r)|| = exp(-D(s')) ^ ||R(r) - R'(r)|| = rexp(-D(s')) = exp(-D(s')) (68) 


for s' > 0{T). Hence, 




*r-|-S 


exp(zi7^*T) = exp iHi^T + i ^ HiT exp i ^ HiT + iH^^T I/'(r) + exp(-D(s')). 


l=il-s 


l=j+l 


(69) 


Using this decomposition sequentially for j = z; —l-|-2ms', m = 0,1, 2,... , mo with mo 
we obtain 


f’r—ll I 1 

2s' 


mo-l-1 mo 

exp(zi7-*r) = n n ^m+ m.oexp(-D(s')), 

m=0 m=0 

i ^ ^ HiT 

l=ii-\-2(m,— l)s' 

( il-2 \ / v+s \ 

+ i Y, HiTU Umo+i = exp z HiT + iH^T , (70) 

l=il—s j y l=ii+2mos' j 

and Vm is an operator acting on = OmiV nr=o Vm with 

s' = 0{T -I-log (mo/C')) such that (i) holds. Clearly, the bond dimension of Kt across the cut j\j + l 
for each z; — 1 < j < /c is upper bounded by 2'^^* ^ = 2'^^^)/poly((5^'). □ 

Lemma 21. There is an MPO K such that (i) 1177 —H|| < z/ed/lOOO =: r]'; (ii) the bond dimension 
across the cuts j\j + 1 for ii — 1 < j < k is upper bounded by D = (1/(5)'^(i+\/*°g(i/’7)/; 

(Hi) ’s in the decomposition for K can be constructed in time poly(D). 

Proof. Following [S], we truncate and discretize the integral in the definition of A. Let tj = rj 

with j = 0, ±1,..., zbT/r, and replace the integral by a sum over tj. The truncation error is 

of order exp for T = 0{^/q log(l/r/')/e'). The discretization error is of or¬ 
der e' t\\H-^\\T/ J2-Kq < n'/S for r = O ( , : V For each tj, Lemma [20] implies that the 

VV^og(im)/ 

propagator exp(—z(i7-* — e'Q)tj) can be approximated to error zy'/S by an MPO of bond dimen¬ 
sion 2'^fo9 poly((5/r/') = 2'^^'^) poly((5/r/') = poly{d/r]') across the cuts j\j -|- 1 

for ii — 1 < j < k. The number of terms is 2T/t + 1 = 0 ^ ^ ^ • Hence the bond 
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dimension across the cut j\j + 1 of the MPO K for ii — 1 < j < k is upper bounded by 
D = (1/(5)*^^^^i°g(i/'5)A), xhe efficiency of constructing K^-^'s in the decomposition 
(j65h for K follows from Lemma [2nlf iii) . □ 

Lemmas dni [2T] imply {(jJ+ rje' for {(j)') = iL|<l)')/||iL|<h')||. It is easy to 
see that Sk '■= ^ is a {k, D^pi, Dp 2 , Ae = r/)-support set. Hence, 

P3 = (i/(5)0(i+Vi°g(i/^)/i°g(i/'5)A) = (l/(5)'^(i). 

We briefly comment on the assumption that we have an estimate Cq of to error Since 

0 < < 1/5 + 2s + 2t + 1 = 0(1/5), we run the whole algorithm (for H-^) with Cq = 

and obtain a “candidate e/’ (cf. Lemma[7]) for each j = 0,1,... , 0{5~^/^). The minimum of these 
candidates is identified as the “final Cj.” This completes our algorithm. 


Acknowledgment 

We would like to thank Xie Chen and Aram W. Harrow for discussions. In particular, we learned 
Lemma [T] from A.W.H some time ago. This work was supported by DARPA OLE. 


References 

[1] D. Aharonov, 1. Arad, and S. Irani. Efficient algorithm for approximating one-dimensional 
ground states. Physical Review A, 82(1):012315, 2010. 

[2] D. Aharonov, 1. Arad, and T. Vidick. Guest column: The quantum PCP conjecture. SIGACT 
News, 44(2):47-79, 2013. 

[3] D. Aharonov, D. Gottesman, S. Irani, and J. Kempe. The power of quantum systems on a 
line. In Proceedings of the fSth Annual IEEE Symposium on Foundations of Computer Science, 
pages 373-383, 2007. 

[4] D. Aharonov, D. Gottesman, S. Irani, and J. Kempe. The power of quantum systems on a 
line. Communications in Mathematical Physics, 287(l):41-65, 2009. 

[5] 1. Arad, A. Kitaev, Z. Landau, and U. Vazirani. An area law and sub-exponential algorithm 
for ID systems. arXiv:1301.1162vl. 

[6] F. G. S. L. Brandao and A. W. Harrow. Product-state approximations to quantum ground 
states. In Proceedings of the 45th Annual ACM Symposium on Theory of Computing, pages 
871-880, 2013. 

[7] C. T. Chubb and S. T. Flammia. Computing the degenerate ground space of gapped spin 
chains in polynomial time. arXiv: 1502.06967. 

[8] M. Fannes, B. Nachtergaele, and R. Werner. Finitely correlated states on quantum spin chains. 
Communications in Mathematical Physics, 144(3) :443-490, 1992. 

[9] M. B. Hastings. Quantum adiabatic computation with a constant gap is not useful in one 
dimension. Physical Review Letters, 103(5) :050502, 2009. 

[10] Y. Huang. Area law in one dimension: Degenerate ground states and Renyi entanglement 
entropy. arXiv: 1403.0327. 


18 




[11] Y. Huang. A polynomial-time algorithm for the ground state of one-dimensional gapped 
Hamiltonians. arXiv; 1406.6355. 

[12] Z. Landau, U. Vazirani, and T. Vidick. A polynomial-time algorithm for the ground state of 
ID gapped local Hamiltonians. arXiv:1307.5143vl. 

[13] T. J. Osborne. Efficient approximation of the dynamics of one-dimensional quantum spin 
systems. Physical Review Letters, 97(15):157202, 2006. 

[14] T. J. Osborne. Hamiltonian complexity. Reports on Progress in Physics, 75(2):022001, 2012. 

[15] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac. Matrix product state representa¬ 
tions. Quantum Information and Computation, 7(5):401~430, 2007. 

[16] N. Schuch and J. I. Cirac. Matrix product state and mean-field solutions for one-dimensional 
systems can be found efficiently. Physical Review A, 82(1);012314, 2010. 

[17] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac. Entropy scaling and simulability by 
matrix product states. Physical Review Letters, 100(3) :030504, 2008. 

[18] F. Verstraete and J. I. Cirac. Matrix product states represent ground states faithfully. Physical 
Review B, 73(9);094423, 2006. 

[19] S. R. White. Density matrix formulation for quantum renormalization groups. Physical Review 
Letters, 69(19):2863-2866, 1992. 

[20] S. R. White. Density-matrix algorithms for quantum renormalization groups. Physical Review 
B, 48(14):10345-10356, 1993. 


19 



